
# ------------------------------------------------------------------------- #
# Replication code for "Covid 19 and Preferences for Progressive Taxation"
# Journal of Elections, Public Opinion and Parties
# Mark Williamson 
# ------------------------------------------------------------------------- #

# PRELIMINARIES ----------------------------------------------

## ---- Load packages ----
library(dplyr); library(estimatr); library(ggplot2);
library(foreign); library(margins); library(stargazer); 
library(xtable); library(lubridate); library(tidyr); library(interflex);
library(sandwich); library(lmtest); library(car); library(sf); library(spdep); 
library(spatialreg)

main.theme <- theme_bw() + 
  theme(legend.key = element_rect(fill = NA, color = NA), 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        text = element_text(family = "Helvetica", size = 18))
theme_set(main.theme)

## ---- Read in data ----

# Load in zip-code level data 
df <- read.csv("./data/main_data/master_zip_data_il.csv")

# Load in ANES data 
anes <- read.csv("./data/main_data/master_anes_data.csv")


# TABLE 1 ------------------------------------------------

## Illinois zip code level --------------------------------

# Define covariate conditioning set
cov <- paste0(c("median_hh_", "gini_index", "pct_latino", 
                "pct_black", "pct_povert", "pct_emp_ed_health", "pct_emp_hospitality",
                "pct_emp_govt", "pct_bachel", "under18", "over65",
                "pop_densit", "ln_pop", "pritz18pct"), 
              collapse = " + ")

# Run models 
modelform <- formula(paste("yestaxpct ~ cases_1k + factor(COUNTYNAME) + ", cov, sep = ""))
mod1 <- lm_robust(modelform, data = df, se_type="HC2")
mod1b <- lm(modelform, data = df)

# Effect of 1 s.d. increase
mod1$coefficients["cases_1k"] * sd(df$cases_1k, na.rm = TRUE)

# Case count increase needed to reach simple majority 
(50-46.7) / mod1$coefficients["cases_1k"]

# Overall statewide cases per 1,000: 
1000 * sum(df$casecount_mail, na.rm = TRUE) / sum(df$total_pop, na.rm = TRUE)

# Percent increase needed 
129.6765 / 21.5

## Chicago deaths --------------------------------

# Define covariate conditioning set
cov <- paste0(c("median_hh_", "gini_index", "pct_latino", 
                "pct_black", "pct_povert", "pct_emp_ed_health", "pct_emp_hospitality",
                "pct_emp_govt", "pct_bachel", "under18", "over65",
                "pop_densit", "ln_pop", "pritz18pct"), 
              collapse = " + ")

modelform <- formula(paste("yestaxpct ~ deaths_1k + ", cov, sep = ""))
mod2 <- lm_robust(modelform, data = df, se_type="HC2")
mod2b <- lm(modelform, data = df)

## Employment and earnings change -----------------------------
modelform <- formula(paste("yestaxpct ~ emp_pct_loss + factor(COUNTYNAME) + ", cov, sep = ""))
mod3 <- lm_robust(modelform, data = df, se_type="HC2")
mod3b <- lm(modelform, data = df)

modelform <- formula(paste("yestaxpct ~ pay_pct_loss + factor(COUNTYNAME) + ", cov, sep = ""))
mod4 <- lm_robust(modelform, data = df, se_type="HC2")
mod4b <- lm(modelform, data = df)

## Lockdowns ---------- 
# (No county FE b/c restrictions are constant within counties)
modelform <- formula(paste("yestaxpct ~ treated_lock + border_seg_fe + ", cov, sep = ""))
mod5 <- lm_robust(modelform, data = df, se_type="HC2")
mod5b <- lm(modelform, data = df)

# Produce table
stargazer::stargazer(mod1b, mod2b, mod3b, mod4b, mod5b,
                     se = c(starprep(mod1), starprep(mod2),
                            starprep(mod3), starprep(mod4),
                            starprep(mod5)), 
                    type = "text", 
                     keep = c("case", "death", "lock", "mean", "emp_pct", "pay_pct"),
                     star.cutoffs = 0.05,
                     omit.stat = "f")


# TABLE 2 ------------------------------------------------------

# Convert panel to first differences
anes_fd <- anes %>% 
  # Those observed in both waves
  filter(panel==1) %>%
  # Rescale financial worry, create egalitarianism index
  mutate(fin_worry_z = (fin_worry - mean(fin_worry[year==2016], na.rm = TRUE)) / 
           sd(fin_worry[year==2016], na.rm = TRUE),
         egal_idx_sd2016 = sd(egal_idx[year==2016], na.rm = TRUE), 
         egal1_sd2016 = sd(egal1[year==2016], na.rm = TRUE), 
         egal2_sd2016 = sd(egal2[year==2016], na.rm = TRUE), 
         egal3_sd2016 = sd(egal3[year==2016], na.rm = TRUE), 
         egal4_sd2016 = sd(egal4[year==2016], na.rm = TRUE) 
  ) %>%
  group_by(id, pid16) %>%
  summarise(mil_tax_fd = mil_tax[year==2020] - mil_tax[year==2016],
            yes_mil_tax_16 = yes_mil_tax[year==2016],
            yes_mil_tax_fd = yes_mil_tax[year==2020] - yes_mil_tax[year==2016],
            egal_idx_fd = (egal_idx[year==2020] - egal_idx[year==2016]),
            egal_idx_fd_z = (egal_idx[year==2020] - egal_idx[year==2016]) / egal_idx_sd2016,
            egal1_fd_z = (egal1[year==2020] - egal1[year==2016]) / egal1_sd2016,
            egal2_fd_z = (egal2[year==2020] - egal2[year==2016]) / egal2_sd2016,
            egal3_fd_z = (egal3[year==2020] - egal3[year==2016]) / egal3_sd2016,
            egal4_fd_z = (egal4[year==2020] - egal4[year==2016]) / egal4_sd2016,
            covid_infect = max(covid_infect),
            chg_working = working[year==2020] - working[year==2016],
            chg_worry =  fin_worry_z[year==2020] - fin_worry_z[year==2016]
  ) %>%
  ungroup() %>%
  distinct()

# Estimate models -- Support for Millionaire's Tax 
mod1 <- lm_robust(yes_mil_tax_fd ~ covid_infect,
                  data=anes_fd, se_type = "HC2")
mod1b <- lm(yes_mil_tax_fd ~ covid_infect,
            data=anes_fd)

mod2 <- lm_robust(yes_mil_tax_fd ~ chg_worry,
                  data=anes_fd, se_type = "HC2")
mod2b <- lm(yes_mil_tax_fd ~ chg_worry,
            data=anes_fd)

mod3 <- lm_robust(yes_mil_tax_fd ~ chg_working,
                  data=anes_fd, se_type = "HC2")
mod3b <- lm(yes_mil_tax_fd ~ chg_working,
            data=anes_fd)

mod4 <- lm_robust(yes_mil_tax_fd ~ covid_infect + chg_worry + chg_working,
                  data=anes_fd, se_type = "HC2")
mod4b <- lm(yes_mil_tax_fd ~ covid_infect + chg_worry + chg_working,
            data=anes_fd)

stargazer::stargazer(mod1b, mod2b, mod3b, mod4b,
                     se = c(starprep(mod1), starprep(mod2),
                            starprep(mod3), starprep(mod4)), 
                     type = "text", 
                     omit.stat = "f")

# APPENDIX -------------------------------------------------------

## FIGURE A1 ------------------------------------------------

### Cases ------------------------------------------------------------
state_pop <- covidcast::state_census %>%
  select(ABBR, POPESTIMATE2019) %>%
  rename(state=ABBR, pop=POPESTIMATE2019)

cdc <- read.csv("./data/sm_data//Weekly_United_States_COVID-19_Cases_and_Deaths_by_State_-_ARCHIVED_20250409.csv") %>%
  mutate(end_date = as.Date(end_date, "%m/%d/%Y")) %>%
  filter(end_date <= as.Date("2020-11-04")) %>%
  left_join(state_pop) %>%
  mutate(tot_cases_1k = tot_cases / pop * 1000,
         tot_deaths_1k = tot_deaths / pop * 1000)

ggplot(cdc %>% filter(end_date==as.Date("2020-11-04"))) + 
  geom_text(aes(tot_cases_1k, tot_deaths_1k, label=state), 
            col="gray70") + 
  geom_text(data=cdc %>% filter(end_date==as.Date("2020-11-04") & 
                                  state=="IL"),
            aes(tot_cases_1k, tot_deaths_1k, label=state), 
            col="black") + 
  labs(x="Pre-election\n Covid-19 cases per\n1,000 residents",
       y="Pre-election\n Covid-19 deaths per\n1,000 residents") + 
  theme(axis.title.y = element_text(angle=0, vjust=0.5))


### Economic impact ------------------------------------------------------------

# GDP change 
bea <- read.csv("./data/sm_data/bea_state_gdp.csv", nrows = 60) %>%
  filter(GeoFips %in% 1:57000 & GeoFips != 11000) %>%
  mutate(q4_q3_chg = (X2020.Q3 - X2019.Q4) / X2019.Q4) %>%
  mutate(GeoName = factor(GeoName, levels = GeoName[order(q4_q3_chg)]))

ggplot(bea) + 
  geom_bar(aes(GeoName, q4_q3_chg), stat="identity", 
           fill="gray70") + 
  geom_bar(data=bea %>% filter(GeoName=="Illinois"),
           aes(GeoName, q4_q3_chg), stat="identity") + 
  theme(axis.text.x = element_text(angle=90, vjust=0.5),
        axis.title.y = element_text(angle=0, vjust=0.5)) + 
  labs(x="", y="Change in\nreal GDP\n from 2019 Q4\nto 2020 Q3") + 
  scale_y_continuous(labels=scales::percent)

# Alternative: Employment change
cty_pop <- covidcast::county_census %>%
  select(FIPS, POPESTIMATE2019, STNAME) %>%
  rename(fips=FIPS, 
         pop=POPESTIMATE2019, 
         state=STNAME) %>%
  mutate(fips = ifelse(substr(fips,1,1)=="0", as.numeric(substr(fips, 2, nchar(fips))), 
                       as.numeric(fips)))

emp_df <- read.csv("./data/sm_data/monthly_emp_chg_by_cty_20_21.csv") %>%
  # Average change in early pandemic vs. over full period 
  group_by(fips, name) %>%
  summarise(emplvl_pct_chg_preelect = mean(emplvl_pct_chg[!month %in% c("October", "November", "December", 
                                                                        "January", "February")], na.rm = TRUE)) %>%
  left_join(cty_pop) %>%
  group_by(state) %>%
  summarise(emplvl_pct_chg_preelect = sum(emplvl_pct_chg_preelect * pop) / sum(pop)) %>%
  mutate(state = factor(state, levels = state[order(emplvl_pct_chg_preelect)]))

ggplot(emp_df) + 
  geom_bar(aes(state, emplvl_pct_chg_preelect/100), stat="identity", 
           fill="gray70") + 
  geom_bar(data=emp_df %>% filter(state=="Illinois"),
           aes(state, emplvl_pct_chg_preelect/100), stat="identity") + 
  theme(axis.text.x = element_text(angle=90, vjust=0.5),
        axis.title.y = element_text(angle=0, vjust=0.5)) + 
  labs(x="", y="Average\nyear-over-year\npre-election\nemployment\nchange") + 
  scale_y_continuous(labels=scales::percent)

# Percentile
1 - ecdf(emp_df$emplvl_pct_chg_preelect)(-9.06)
rank(emp_df$emplvl_pct_chg_preelect)[which(emp_df$state=="Illinois")]


## FIGURE 2 -----------------------------------------------
anes20 <- read.csv("./data/sm_data/anes_timeseries_2020_csv_20220210.csv") %>%
  mutate(dem_therm = V201156,
         rep_therm = V201157, 
         pid = case_when(V201231x %in% 1:3 ~ "Democrat", 
                         V201231x %in% 5:7 ~ "Republican"),
         in_therm = ifelse(pid=="Democrat", dem_therm, rep_therm),
         out_therm = ifelse(pid=="Democrat", rep_therm, dem_therm),
         aff_pol = in_therm - out_therm,
         state = case_when(
           V203000 == 1  ~ "Alabama",
           V203000 == 2  ~ "Alaska",
           V203000 == 4  ~ "Arizona",
           V203000 == 5  ~ "Arkansas",
           V203000 == 6  ~ "California",
           V203000 == 8  ~ "Colorado",
           V203000 == 9  ~ "Connecticut",
           V203000 == 10 ~ "Delaware",
           V203000 == 11 ~ NA,
           V203000 == 12 ~ "Florida",
           V203000 == 13 ~ "Georgia",
           V203000 == 15 ~ "Hawaii",
           V203000 == 16 ~ "Idaho",
           V203000 == 17 ~ "Illinois",
           V203000 == 18 ~ "Indiana",
           V203000 == 19 ~ "Iowa",
           V203000 == 20 ~ "Kansas",
           V203000 == 21 ~ "Kentucky",
           V203000 == 22 ~ "Louisiana",
           V203000 == 23 ~ "Maine",
           V203000 == 24 ~ "Maryland",
           V203000 == 25 ~ "Massachusetts",
           V203000 == 26 ~ "Michigan",
           V203000 == 27 ~ "Minnesota",
           V203000 == 28 ~ "Mississippi",
           V203000 == 29 ~ "Missouri",
           V203000 == 30 ~ "Montana",
           V203000 == 31 ~ "Nebraska",
           V203000 == 32 ~ "Nevada",
           V203000 == 33 ~ "New Hampshire",
           V203000 == 34 ~ "New Jersey",
           V203000 == 35 ~ "New Mexico",
           V203000 == 36 ~ "New York",
           V203000 == 37 ~ "North Carolina",
           V203000 == 38 ~ "North Dakota",
           V203000 == 39 ~ "Ohio",
           V203000 == 40 ~ "Oklahoma",
           V203000 == 41 ~ "Oregon",
           V203000 == 42 ~ "Pennsylvania",
           V203000 == 44 ~ "Rhode Island",
           V203000 == 45 ~ "South Carolina",
           V203000 == 46 ~ "South Dakota",
           V203000 == 47 ~ "Tennessee",
           V203000 == 48 ~ "Texas",
           V203000 == 49 ~ "Utah",
           V203000 == 50 ~ "Vermont",
           V203000 == 51 ~ "Virginia",
           V203000 == 53 ~ "Washington",
           V203000 == 54 ~ "West Virginia",
           V203000 == 55 ~ "Wisconsin",
           V203000 == 56 ~ "Wyoming",
           TRUE ~ NA_character_  # Handles unexpected codes
         )
  ) %>%
  filter(!is.na(pid)) %>%
  select(state, pid, dem_therm, rep_therm, aff_pol) %>%
  group_by(state) %>%
  summarise(n=n(), 
            aff_pol_avg = mean(aff_pol, na.rm = TRUE),
            aff_pol_se = sd(aff_pol, na.rm = TRUE) / sqrt(n)) %>%
  filter(!is.na(state)) %>%
  mutate(state = factor(state, levels = state[order(aff_pol_avg)]))

# Percentile
ecdf(anes20$aff_pol_avg)(48.18710)
rank(anes20$aff_pol_avg)[which(anes20$state=="Illinois")]


ggplot(anes20) + 
  geom_point(aes(state, aff_pol_avg), col="gray70") + 
  geom_linerange(aes(state, ymin = aff_pol_avg - 1.96*aff_pol_se,
                     ymax = aff_pol_avg + 1.96*aff_pol_se),
                 col="gray70") + 
  geom_point(data=anes20 %>% filter(state %in% c("Illinois")),
             aes(state, aff_pol_avg), col="black") + 
  geom_linerange(data=anes20 %>% filter(state %in% c("Illinois")),
                 aes(state, ymin = aff_pol_avg - 1.96*aff_pol_se,
                     ymax = aff_pol_avg + 1.96*aff_pol_se),
                 col="black") + 
  theme(axis.text.x = element_text(angle=90, vjust=0.5),
        axis.title.y = element_text(angle=0, vjust=0.5)) + 
  labs(x="", y="Affective\nPolarization")


## TABLE A1 -----------------------------------------

# Create a new variable for cases/1k before E.D. 
df <- df %>% 
  mutate(cases_1k_ed = (casecount_ed / total_pop)*1000)

# Define covariate conditioning set
cov <- paste0(c("median_hh_", "gini_index", "pct_latino", 
                "pct_black", "pct_povert", "pct_emp_ed_health", "pct_emp_hospitality",
                "pct_emp_govt", "pct_bachel", "under18", "over65",
                "pop_densit", "ln_pop", "pritz18pct"), 
              collapse = " + ")

# Run models 
modelform <- formula(paste("yestaxpct ~ cases_1k + factor(COUNTYNAME) + ", cov, sep = ""))
mod1 <- lm_robust(modelform, data = df, se_type="HC2")
mod1b <- lm(modelform, data = df)

modelform <- formula(paste("yestaxpct ~ cases_1k_ed + factor(COUNTYNAME) + ", cov, sep = ""))
mod2 <- lm_robust(modelform, data = df, se_type="HC2")
mod2b <- lm(modelform, data = df)

stargazer::stargazer(mod1b, mod2b, 
                     se = c(starprep(mod1), starprep(mod2)), 
                        type = "text", 
                     keep = c("case", "death", "lock", "mean"),
                     star.cutoffs = 0.05,
                     omit.stat = "f")


## FIGURE A3 -----------------------------------------
emp_df <- read.csv("./data/sm_data/monthly_emp_chg_by_cty_20_21.csv") %>%
  # Average change in early pandemic vs. over full period 
  group_by(fips, name) %>%
  summarise(emplvl_pct_chg_preelect = mean(emplvl_pct_chg[!month %in% c("October", "November", "December", 
                                                                        "January", "February")], na.rm = TRUE),
            emplvl_pct_chg_postelect = mean(emplvl_pct_chg[month %in% c("October", "November", "December", 
                                                                        "January", "February")], na.rm = TRUE)) %>%
  filter(grepl("Illinois", name))

cor(emp_df$emplvl_pct_chg_postelect, emp_df$emplvl_pct_chg_preelect)

ggplot(emp_df) + 
  geom_point(aes(emplvl_pct_chg_preelect, emplvl_pct_chg_postelect), alpha=0.4) + 
  geom_smooth(aes(emplvl_pct_chg_preelect, emplvl_pct_chg_postelect),
              method="lm", col="red", se=F) + 
#  geom_smooth(data=emp_df %>% filter(emplvl_pct_chg_preelect > -20), 
#              aes(emplvl_pct_chg_preelect, emplvl_pct_chg_postelect),
#              method="lm", col="blue", se=F) + 
  geom_text(data=data.frame(), aes(-25, -5), label="r=0.72", size=5) + 
  labs(x="Average year-over-year\nemployment change\nbefore 2020 Election",
       y="Average\nyear-over-year\nemployment\nchange after\n2020 Election") + 
  theme(axis.title.y = element_text(angle=0, vjust=0.5)) + 
  scale_x_continuous(breaks=seq(-30,0,5), limits=c(-30,0))

## TABLE A2 --------------------------

plot_df <- df %>% 
  select(c("median_hh_", "gini_index", "pct_latino", 
           "pct_black", "pct_povert", "pct_emp_ed_health", "pct_emp_hospitality",
           "pct_emp_govt", "pct_bachel", "under18", "over65",
           "pop_densit", "total_pop", "pritz18pct", "cases_1k"), 
         treated_lock) %>%
  filter(!is.na(treated_lock)) %>%
  group_by(treated_lock) %>% 
  summarise_all(.funs=mean) %>%
  t()

print(xtable(plot_df, type = "latex"))
t.test(df$cases_1k[df$treated_lock==1], df$cases_1k[df$treated_lock==0])

## TABLE A3 ----------------------------------

# No covariates
mod1 <- lm_robust(yestaxpct ~ treated_lock + border_seg_fe, data = df, se_type="HC2")
mod1b <- lm(yestaxpct ~ treated_lock + border_seg_fe, data = df)

# Covariates 
modelform <- formula(paste("yestaxpct ~ treated_lock + border_seg_fe + ", cov, sep = ""))
mod2 <- lm_robust(modelform, data = df, se_type="HC2")
mod2b <- lm(modelform, data = df)

# Removing zip codes that straddle the border
modelform <- formula(paste("yestaxpct ~ treated_lock + border_seg_fe + ", cov, sep = ""))
mod3 <- lm_robust(modelform, data = df %>% filter(closecall_lock==0), se_type="HC2")
mod3b <- lm(modelform, data = df %>% filter(closecall_lock==0))

# Employment loss
modelform <- formula(paste("emp_pct_loss ~ treated_lock + border_seg_fe + ", cov, sep = ""))
mod4 <- lm_robust(modelform, data = df, se_type="HC2")
mod4b <- lm(modelform, data = df)

# Payroll loss 
modelform <- formula(paste("pay_pct_loss ~ treated_lock + border_seg_fe + ", cov, sep = ""))
mod5 <- lm_robust(modelform, data = df, se_type="HC2")
mod5b <- lm(modelform, data = df)

# Produce table
stargazer::stargazer(mod4b, mod5b, mod1b, mod2b, mod3b, 
                     se = c(starprep(mod4), starprep(mod5), starprep(mod1),
                            starprep(mod2), starprep(mod3)), 
                          type = "text", 
                     keep = c("case", "death", "lock", "mean"),
                     star.cutoffs = 0.05,
                     omit.stat = "f")


## FIGURE A6 ---------------------------------

#### Case counts --------------------------------

# Restrict to available observations
cate_df <- df %>% filter(!is.na(pritz18pct)) 

out <- interflex(estimator = "binning",
                 vcov.type = "robust",
                 Y="yestaxpct", D="cases_1k", X = "pritz18pct",
                 Z=c("median_hh_", "gini_index", "pct_latino", 
                     "pct_black", "pct_povert", "pct_emp_ed_health", "pct_emp_hospitality",
                     "pct_emp_govt", "pct_bachel", "under18", "over65",
                     "pop_densit", "ln_pop"), full.moderate = TRUE,
                 FE = "COUNTYNAME",
                 data=cate_df, 
                 na.rm = TRUE)

ggplot() + 
  geom_hline(yintercept = 0, lty="dashed", col="grey") + 
  geom_rug(aes(cate_df$pritz18pct), alpha=0.3) + 
  geom_line(data=out$est.lin$`cases_1k=13.458 (50%)`, 
            aes(X, ME)) + 
  geom_ribbon(data=out$est.lin$`cases_1k=13.458 (50%)`, 
              aes(X, ymin=`lower CI(95%)`, ymax=`upper CI(95%)`), 
              alpha=0.3) + 
  geom_pointrange(data=out$est.bin$`cases_1k=13.458 (50%)`, 
                  aes(x=x0, y=coef, ymin=CI.lower, ymax=CI.upper), stat="identity") + 
  scale_x_continuous(limits = c(0,1)) + 
  labs(x="\nDemocratic gubernatorial\nvote share (2018)", 
       y="Implied\neffect of\nCovid-19 cases\nper 1,000 on\nsupport for\nprogressive\ntaxation") + 
  theme(axis.title.y = element_text(angle = 0, vjust=0.5)) 

#### Employment loss ----------------------------

# Restrict to available observations
cate_df <- df %>% filter(!is.na(pritz18pct)) 

out <- interflex(estimator = "binning",
                 vcov.type = "robust",
                 Y="yestaxpct", D="emp_pct_loss", X = "pritz18pct",
                 Z=c("median_hh_", "gini_index", "pct_latino", 
                     "pct_black", "pct_povert", "pct_emp_ed_health", "pct_emp_hospitality",
                     "pct_emp_govt", "pct_bachel", "under18", "over65",
                     "pop_densit", "ln_pop"), full.moderate = TRUE,
                 FE = "COUNTYNAME",
                 data=cate_df, 
                 na.rm = TRUE)

ggplot() + 
  geom_hline(yintercept = 0, lty="dashed", col="grey") + 
  geom_rug(aes(cate_df$pritz18pct), alpha=0.3) + 
  geom_line(data=out$est.lin$`emp_pct_loss=3.412 (50%)`, 
            aes(X, ME)) + 
  geom_ribbon(data=out$est.lin$`emp_pct_loss=3.412 (50%)`, 
              aes(X, ymin=`lower CI(95%)`, ymax=`upper CI(95%)`), 
              alpha=0.3) + 
  geom_pointrange(data=out$est.bin$`emp_pct_loss=3.412 (50%)`, 
                  aes(x=x0, y=coef, ymin=CI.lower, ymax=CI.upper), stat="identity") + 
  scale_x_continuous(limits = c(0,1)) + 
  labs(x="\nDemocratic gubernatorial\nvote share (2018)", 
       y="Implied effect\npercent of\nemployment\nlosses on\nsupport for\nprogressive\ntaxation") + 
  theme(axis.title.y = element_text(angle = 0, vjust=0.5)) + 
  guides(col="none")


#### Earnings loss ----------------------------

# Restrict to available observations
cate_df <- df %>% filter(!is.na(pritz18pct)) 

out <- interflex(estimator = "binning",
                 vcov.type = "robust",
                 Y="yestaxpct", D="pay_pct_loss", X = "pritz18pct",
                 Z=c("median_hh_", "gini_index", "pct_latino", 
                     "pct_black", "pct_povert", "pct_emp_ed_health", "pct_emp_hospitality",
                     "pct_emp_govt", "pct_bachel", "under18", "over65",
                     "pop_densit", "ln_pop"), full.moderate = TRUE,
                 FE = "COUNTYNAME",
                 data=cate_df, 
                 na.rm = TRUE)

ggplot() + 
  geom_hline(yintercept = 0, lty="dashed", col="grey") + 
  geom_rug(aes(cate_df$pritz18pct), alpha=0.3) + 
  geom_line(data=out$est.lin$`pay_pct_loss=-6.005 (50%)`, 
            aes(X, ME)) + 
  geom_ribbon(data=out$est.lin$`pay_pct_loss=-6.005 (50%)`, 
              aes(X, ymin=`lower CI(95%)`, ymax=`upper CI(95%)`), 
              alpha=0.3) + 
  geom_pointrange(data=out$est.bin$`pay_pct_loss=-6.005 (50%)`, 
                  aes(x=x0, y=coef, ymin=CI.lower, ymax=CI.upper), stat="identity") + 
  scale_x_continuous(limits = c(0,1)) + 
  labs(x="\nDemocratic gubernatorial\nvote share (2018)", 
       y="Implied effect\npercent of\nemployment\nlosses on\nsupport for\nprogressive\ntaxation") + 
  theme(axis.title.y = element_text(angle = 0, vjust=0.5)) + 
  guides(col="none")


## FIGURE A7 ---------------------------------
(mod1 <- lm_robust(yes_mil_tax_fd ~ covid_infect:pid16,
                   data=anes_fd, se_type = "HC2"))
(mod2 <- lm_robust(yes_mil_tax_fd ~ chg_worry:pid16,
                   data=anes_fd, se_type = "HC2"))
(mod3 <- lm_robust(yes_mil_tax_fd ~ chg_working:pid16,
                   data=anes_fd, se_type = "HC2"))

margfx <- rbind(margins(mod1, variables = "covid_infect", 
                        at = list(pid16=unique(anes$pid16[!is.na(anes$pid16)]))) %>% 
                  summary() %>% 
                  mutate(var="Covid-19\ninfection"),
                margins(mod2, variables = "chg_worry", 
                        at = list(pid16=unique(anes$pid16[!is.na(anes$pid16)]))) %>% 
                  summary() %>% 
                  mutate(var="Financial\nworry"),
                margins(mod3, variables = "chg_working", 
                        at = list(pid16=unique(anes$pid16[!is.na(anes$pid16)]))) %>% 
                  summary() %>% 
                  mutate(var="Working\nstatus")) %>%
  mutate(pid16 = factor(pid16, levels=c("Democrat", "Independent", "Republican"),
                        ordered = TRUE))

ggplot(margfx) + 
  facet_wrap(~var) + 
  geom_hline(yintercept = 0, lty="dashed", col="grey") + 
  geom_point(aes(pid16, AME, col=pid16)) + 
  geom_linerange(aes(pid16, ymin=lower, max=upper, col=pid16)) + 
  scale_color_manual(values=c("blue", "gray30", "red")) + 
  labs(x="", y="Implied effect\non support for\nprogressive\ntaxation\n(0/1)") + 
  guides(col="none") + 
  theme(axis.title.y = element_text(angle=0, vjust=0.5)) 

## FIGURE A8 ------------------------------------------

# Set up list of all covariates used in main model (except county FE) 
cov <- c("median_hh_", "gini_index", "pct_latino", 
         "pct_black", "pct_povert", "pct_emp_ed_health", "pct_emp_hospitality",
         "pct_emp_govt", "pct_bachel", "under18", "over65",
         "pop_densit", "ln_pop", "pritz18pct")
cov_lab <- c("Median\nHH income", "Gini\nindex", "Latino %", "Black %",
             "% in\npoverty", "% in Ed.,\nhealth", "% in\nhospitality",
             "% in\ngovt.", "% univsty", "% under\n18", "% over\n65",
             "Pop.\ndensity", "Log pop.", "Dem.\nvote share")

### Case counts ----------------------------------------------
# Set up vector to store model estimates
est <- vector()
se <- vector() 

for (i in 1:length(cov)) {
  
  # Determine relevant covariates 
  cov_temp <- cov[1:length(cov) != i]
  
  # Paste together new covariate list 
  cov_temp <- paste0(cov_temp, collapse = " + ")
  
  # Run model using all but the excluded covariate 
  modelform <- formula(paste("yestaxpct ~ cases_1k + factor(COUNTYNAME) + ", 
                             cov_temp, sep = ""))
  mod <- lm_robust(modelform, data = df, se_type="HC2")
  
  est[i] <- mod$coefficients["cases_1k"]
  se[i] <- mod$std.error["cases_1k"]
}

# plot estimates 
ggplot() + 
  geom_hline(yintercept = 0, lty="dotted", col="red") + 
  geom_point(aes(1:length(est), est)) + 
  geom_segment(aes(x = 1:length(est), xend = 1:length(est),
                   y=est-1.96*se, yend=est+1.96*se)) + 
  scale_x_continuous(labels = cov_lab, breaks=1:length(est)) + 
  labs(x="\nExcluded\ncovariate", y = "Coefficient\nestimate \n on Covid cases\nper 1,000") + 
  theme(axis.text.x = element_text(angle = 0),
        axis.title.y = element_text(angle = 0, vjust = 0.5))

### Employment losses ----------------------------------------------
# Set up vector to store model estimates
est <- vector()
se <- vector() 

for (i in 1:length(cov)) {
  
  # Determine relevant covariates 
  cov_temp <- cov[1:length(cov) != i]
  
  # Paste together new covariate list 
  cov_temp <- paste0(cov_temp, collapse = " + ")
  
  # Run model using all but the excluded covariate 
  modelform <- formula(paste("yestaxpct ~ emp_pct_loss + factor(COUNTYNAME) + ", 
                             cov_temp, sep = ""))
  mod <- lm_robust(modelform, data = df, se_type="HC2")
  
  est[i] <- mod$coefficients["emp_pct_loss"]
  se[i] <- mod$std.error["emp_pct_loss"]
}

# plot estimates 
ggplot() + 
  geom_hline(yintercept = 0, lty="dotted", col="red") + 
  geom_point(aes(1:length(est), est)) + 
  geom_segment(aes(x = 1:length(est), xend = 1:length(est),
                   y=est-1.96*se, yend=est+1.96*se)) + 
  scale_x_continuous(labels = cov_lab, breaks=1:length(est)) + 
  labs(x="\nExcluded\ncovariate", 
       y = "Coefficient\nestimate \n on percent\nemployment\nloss") + 
  theme(axis.text.x = element_text(angle = 0),
        axis.title.y = element_text(angle = 0, vjust = 0.5))

### Earnings losses ----------------------------------------------
# Set up vector to store model estimates
est <- vector()
se <- vector() 

for (i in 1:length(cov)) {
  
  # Determine relevant covariates 
  cov_temp <- cov[1:length(cov) != i]
  
  # Paste together new covariate list 
  cov_temp <- paste0(cov_temp, collapse = " + ")
  
  # Run model using all but the excluded covariate 
  modelform <- formula(paste("yestaxpct ~ pay_pct_loss + factor(COUNTYNAME) + ", 
                             cov_temp, sep = ""))
  mod <- lm_robust(modelform, data = df, se_type="HC2")
  
  est[i] <- mod$coefficients["pay_pct_loss"]
  se[i] <- mod$std.error["pay_pct_loss"]
}

# plot estimates 
ggplot() + 
  geom_hline(yintercept = 0, lty="dotted", col="red") + 
  geom_point(aes(1:length(est), est)) + 
  geom_segment(aes(x = 1:length(est), xend = 1:length(est),
                   y=est-1.96*se, yend=est+1.96*se)) + 
  scale_x_continuous(labels = cov_lab, breaks=1:length(est)) + 
  labs(x="\nExcluded\ncovariate", 
       y = "Coefficient\nestimate \n on percent\nearnings\nloss") + 
  theme(axis.text.x = element_text(angle = 0),
        axis.title.y = element_text(angle = 0, vjust = 0.5))


## FIGURE A9 ------------------------------

# Define essential covariate conditioning set 
cov <- paste0(c("median_hh_", "gini_index", "pct_latino", 
                "pct_black", "pct_povert", "pct_emp_ed_health", "pct_emp_hospitality",
                "pct_emp_govt", "pct_bachel", "under18", "over65",
                "pop_densit", "ln_pop", "pritz18pct"), 
              collapse = " + ")

# Set list of counties 
counties <- sort(unique(df$COUNTYNAME))

### Case counts -----------------------------------------

# Set up vector to store model estimates
est <- vector()
se <- vector() 

# Loop over counties 
for (i in 1:length(counties)) {
  
  # Run model using all but the excluded county 
  modelform <- formula(paste("yestaxpct ~ cases_1k + factor(COUNTYNAME) + ", cov, sep = ""))
  mod <- lm_robust(modelform, data = filter(df, COUNTYNAME!=counties[i]), se_type="HC2")
  
  est[i] <- mod$coefficients["cases_1k"]
  se[i] <- mod$std.error["cases_1k"]
}


# plot estimates 
ggplot() + 
  geom_hline(yintercept = 0, lty="dotted", col="red") + 
  geom_point(aes(1:length(est), est)) + 
  geom_segment(aes(x = 1:length(est), xend = 1:length(est),
                   y=est-1.96*se, yend=est+1.96*se)) + 
  scale_x_continuous(labels = counties, breaks=1:length(counties)) + 
  labs(x="Excluded county", y = "Coefficient\nestimate \n on Covid cases\nper 1,000") + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        axis.title.y = element_text(angle=0, vjust=0.5))


### Employment losses -----------------------------------------

# Set up vector to store model estimates
est <- vector()
se <- vector() 

# Loop over counties 
for (i in 1:length(counties)) {
  
  # Run model using all but the excluded county 
  modelform <- formula(paste("yestaxpct ~ emp_pct_loss + factor(COUNTYNAME) + ", cov, sep = ""))
  mod <- lm_robust(modelform, data = filter(df, COUNTYNAME!=counties[i]), se_type="HC2")
  
  est[i] <- mod$coefficients["emp_pct_loss"]
  se[i] <- mod$std.error["emp_pct_loss"]
}


# plot estimates 
ggplot() + 
  geom_hline(yintercept = 0, lty="dotted", col="red") + 
  geom_point(aes(1:length(est), est)) + 
  geom_segment(aes(x = 1:length(est), xend = 1:length(est),
                   y=est-1.96*se, yend=est+1.96*se)) + 
  scale_x_continuous(labels = counties, breaks=1:length(counties)) + 
  labs(x="Excluded county", y = "Coefficient\nestimate \n on percent\nemployment\nloss") + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        axis.title.y = element_text(angle=0, vjust=0.5))


### Earnings losses -----------------------------------------

# Set up vector to store model estimates
est <- vector()
se <- vector() 

# Loop over counties 
for (i in 1:length(counties)) {
  
  # Run model using all but the excluded county 
  modelform <- formula(paste("yestaxpct ~ pay_pct_loss + factor(COUNTYNAME) + ", cov, sep = ""))
  mod <- lm_robust(modelform, data = filter(df, COUNTYNAME!=counties[i]), se_type="HC2")
  
  est[i] <- mod$coefficients["pay_pct_loss"]
  se[i] <- mod$std.error["pay_pct_loss"]
}


# plot estimates 
ggplot() + 
  geom_hline(yintercept = 0, lty="dotted", col="red") + 
  geom_point(aes(1:length(est), est)) + 
  geom_segment(aes(x = 1:length(est), xend = 1:length(est),
                   y=est-1.96*se, yend=est+1.96*se)) + 
  scale_x_continuous(labels = counties, breaks=1:length(counties)) + 
  labs(x="Excluded county", y = "Coefficient\nestimate \n on percent\nearnings\nloss") + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        axis.title.y = element_text(angle=0, vjust=0.5))


## TABLE A4 -------------------------------

# Load in county-level data 
cty <- read.csv("./data/sm_data/master_county_data.csv")

cov <- paste0(c("median_hh_", "gini_index", "pct_latino", 
                "pct_black", "pct_povert", "pct_emp_ed_health", "pct_emp_hospitality",
                "pct_emp_govt", "pct_bachel", "under18", "over65",
                "pop_densit", "ln_pop", "pritz18pct", "region"), 
              collapse = " + ")

# Run models 
modelform <- formula(paste("yestax ~ cases_1k + ", cov, sep = ""))
mod1 <- lm_robust(modelform, data = cty, se_type="HC2")
mod1b <- lm(modelform, data = cty)

modelform <- formula(paste("yestax ~ deaths_1k + ", cov, sep = ""))
mod2 <- lm_robust(modelform, data = cty, se_type="HC2")
mod2b <- lm(modelform, data = cty)

modelform <- formula(paste("yestax ~ mean_yoy_mthly_emplvl_pct_chg + ", cov, sep = ""))
mod3 <- lm_robust(modelform, data = cty, se_type="HC2")
mod3b <- lm(modelform, data = cty)

stargazer::stargazer(mod1b, mod2b, mod3b, 
                     se = c(starprep(mod1), starprep(mod2),
                            starprep(mod3)), 
                     type = "text",
                     star.cutoffs = 0.05,
                     keep = c("case", "death", "days", "mean"),
                     omit.stat = "f")

## TABLE A5 ----------------------------------

# Spatial data
sp_df <- st_read("./data/sm_data/spatial_data/il_combined_zip.shp") %>%
  select(zipcode) %>%
  mutate(zipcode = as.numeric(zipcode)) %>%
  left_join(df, by="zipcode")

# Convert to projection with meters as units 
sp_df <- st_transform(sp_df, crs = 32616)

### Case counts ------------------

# Restrict to non-missing obs
est_df <- sp_df %>%
  filter(if_all(c(
    "median_hh_", "gini_index", "pct_latino", 
    "pct_black", "pct_povert", "pct_emp_ed_health", 
    "pct_emp_hospitality", "pct_emp_govt", "pct_bachel", 
    "under18", "over65", "pop_densit", "total_pop",
    "ln_pop", "pritz18pct", "yestaxpct", "cases_1k"
  ), ~ !is.na(.)))

# All neighbours within d2 km
nb <- dnearneigh(x = st_centroid(est_df), d1 = 0, d2 = 25000)

# Given a list of neighbors, we can use nbdists() to compute the distances along the links. 
# Then, we can construct the list with spatial weights based on inverse distance values using 
# nb2listw() where the argument glist is equal to a list of general weights corresponding to neighbors.
dists <- nbdists(nb, st_centroid(est_df)) 
ids <- lapply(dists, function(x){1/x}) # inverse distances 
nbw <- nb2listw(nb, glist = ids, style = "W", zero.policy=TRUE) # Even though you’re using style = "B" (binary), when you supply a glist of weights, those custom weights override the binary 1s/0s. The "B" style just tells nb2listw() not to row-standardize your custom weights — i.e., your inverse-distance weights are used as-is.

# Moran's I 
moran.test(est_df$cases_1k, nbw, alternative="greater")

# Spatial model can't accomdate county FE, so we will demean before estimation 
vars_to_demean <- c(
  "median_hh_", "gini_index", "pct_latino", "pct_black", "pct_povert",
  "pct_emp_ed_health", "pct_emp_hospitality", "pct_emp_govt",
  "pct_bachel", "under18", "over65", "pop_densit", "ln_pop", "pritz18pct", 
  "yestaxpct", "cases_1k"
)
est_df <- est_df %>%
  group_by(COUNTYNAME) %>%
  mutate(across(
    all_of(vars_to_demean),
    ~ . - mean(., na.rm = TRUE),
    .names = "{.col}_dm"
  )) %>%
  ungroup()

# Estimate using de-meaned Y and Xs
vars_to_est <- paste0(paste0(c(
  "median_hh_", "gini_index", "pct_latino", "pct_black", "pct_povert",
  "pct_emp_ed_health", "pct_emp_hospitality", "pct_emp_govt",
  "pct_bachel", "under18", "over65", "pop_densit", "ln_pop", "pritz18pct"), "_dm"), collapse = " + ")

modelform <- formula(paste("yestaxpct_dm ~ cases_1k_dm + ", vars_to_est, sep = ""))
mod <- lmSLX(modelform, 
             data=est_df, listw=nbw)
summary(impacts(mod)) # Assumes homoskedasticity 

# Calculate robust standard errors
robust_vcov <- vcovHC(mod, type = "HC2")
coeftest(mod, vcov. = robust_vcov)

# Delta method to get robust SE for sum of the two coefficients (total "effect")
(delta_out <- deltaMethod(object = mod$coefficients,
                          g = "cases_1k_dm + lag.cases_1k_dm",
                          vcov. = robust_vcov))


### Employment loss --------------------------------

# Restrict to non-missing obs
est_df <- sp_df %>%
  filter(if_all(c(
    "median_hh_", "gini_index", "pct_latino", 
    "pct_black", "pct_povert", "pct_emp_ed_health", 
    "pct_emp_hospitality", "pct_emp_govt", "pct_bachel", 
    "under18", "over65", "pop_densit", "total_pop",
    "ln_pop", "pritz18pct", "yestaxpct", "emp_pct_loss"
  ), ~ !is.na(.)))

# All neighbours within d2 km
nb <- dnearneigh(x = st_centroid(est_df), d1 = 0, d2 = 25000)

# Given a list of neighbors, we can use nbdists() to compute the distances along the links. 
# Then, we can construct the list with spatial weights based on inverse distance values using 
# nb2listw() where the argument glist is equal to a list of general weights corresponding to neighbors.
dists <- nbdists(nb, st_centroid(est_df)) 
ids <- lapply(dists, function(x){1/x}) # inverse distances 
nbw <- nb2listw(nb, glist = ids, style = "W", zero.policy=TRUE) # Even though you’re using style = "B" (binary), when you supply a glist of weights, those custom weights override the binary 1s/0s. The "B" style just tells nb2listw() not to row-standardize your custom weights — i.e., your inverse-distance weights are used as-is.

# Moran's I 
moran.test(est_df$emp_pct_loss, nbw, alternative="greater")

# Spatial model can't accomdate county FE, so we will demean before estimation 
vars_to_demean <- c(
  "median_hh_", "gini_index", "pct_latino", "pct_black", "pct_povert",
  "pct_emp_ed_health", "pct_emp_hospitality", "pct_emp_govt",
  "pct_bachel", "under18", "over65", "pop_densit", "ln_pop", "pritz18pct", 
  "yestaxpct", "emp_pct_loss"
)
est_df <- est_df %>%
  group_by(COUNTYNAME) %>%
  mutate(across(
    all_of(vars_to_demean),
    ~ . - mean(., na.rm = TRUE),
    .names = "{.col}_dm"
  )) %>%
  ungroup()

# Estimate using de-meaned Y and Xs
vars_to_est <- paste0(paste0(c(
  "median_hh_", "gini_index", "pct_latino", "pct_black", "pct_povert",
  "pct_emp_ed_health", "pct_emp_hospitality", "pct_emp_govt",
  "pct_bachel", "under18", "over65", "pop_densit", "ln_pop", "pritz18pct"), "_dm"), collapse = " + ")


modelform <- formula(paste("yestaxpct_dm ~ emp_pct_loss_dm + ", vars_to_est, sep = ""))
mod <- lmSLX(modelform, 
             data=est_df, listw=nbw)
summary(impacts(mod)) # Assumes homoskedasticity 

# Calculate robust standard errors
robust_vcov <- vcovHC(mod, type = "HC2")
coeftest(mod, vcov. = robust_vcov)

# Delta method to get robust SE for sum of the two coefficients (total "effect")
(delta_out <- deltaMethod(object = mod$coefficients,
                          g = "emp_pct_loss_dm + lag.emp_pct_loss_dm",
                          vcov. = robust_vcov))

### Earnings loss --------------------------------
# Restrict to non-missing obs
est_df <- sp_df %>%
  filter(if_all(c(
    "median_hh_", "gini_index", "pct_latino", 
    "pct_black", "pct_povert", "pct_emp_ed_health", 
    "pct_emp_hospitality", "pct_emp_govt", "pct_bachel", 
    "under18", "over65", "pop_densit", "total_pop",
    "ln_pop", "pritz18pct", "yestaxpct", "pay_pct_loss"
  ), ~ !is.na(.)))

# All neighbours within d2 km
nb <- dnearneigh(x = st_centroid(est_df), d1 = 0, d2 = 25000)

# Given a list of neighbors, we can use nbdists() to compute the distances along the links. 
# Then, we can construct the list with spatial weights based on inverse distance values using 
# nb2listw() where the argument glist is equal to a list of general weights corresponding to neighbors.
dists <- nbdists(nb, st_centroid(est_df)) 
ids <- lapply(dists, function(x){1/x}) # inverse distances 
nbw <- nb2listw(nb, glist = ids, style = "W", zero.policy=TRUE) # Even though you’re using style = "B" (binary), when you supply a glist of weights, those custom weights override the binary 1s/0s. The "B" style just tells nb2listw() not to row-standardize your custom weights — i.e., your inverse-distance weights are used as-is.

# Moran's I 
moran.test(est_df$pay_pct_loss, nbw, alternative="greater")

# Spatial model can't accomdate county FE, so we will demean before estimation 
vars_to_demean <- c(
  "median_hh_", "gini_index", "pct_latino", "pct_black", "pct_povert",
  "pct_emp_ed_health", "pct_emp_hospitality", "pct_emp_govt",
  "pct_bachel", "under18", "over65", "pop_densit", "ln_pop", "pritz18pct", 
  "yestaxpct", "pay_pct_loss"
)
est_df <- est_df %>%
  group_by(COUNTYNAME) %>%
  mutate(across(
    all_of(vars_to_demean),
    ~ . - mean(., na.rm = TRUE),
    .names = "{.col}_dm"
  )) %>%
  ungroup()

# Estimate using de-meaned Y and Xs
vars_to_est <- paste0(paste0(c(
  "median_hh_", "gini_index", "pct_latino", "pct_black", "pct_povert",
  "pct_emp_ed_health", "pct_emp_hospitality", "pct_emp_govt",
  "pct_bachel", "under18", "over65", "pop_densit", "ln_pop", "pritz18pct"), "_dm"), collapse = " + ")

modelform <- formula(paste("yestaxpct_dm ~ pay_pct_loss_dm + ", vars_to_est, sep = ""))
mod <- lmSLX(modelform, 
             data=est_df, listw=nbw)
summary(impacts(mod)) # Assumes homoskedasticity 

# Calculate robust standard errors
robust_vcov <- vcovHC(mod, type = "HC2")
round(coeftest(mod, vcov. = robust_vcov), 3)

# Delta method to get robust SE for sum of the two coefficients (total "effect")

(delta_out <- deltaMethod(object = mod$coefficients,
                          g = "pay_pct_loss_dm + lag.pay_pct_loss_dm",
                          vcov. = robust_vcov))


## TABLE A6 ---------------------------------

# Define covariate conditioning set
cov <- paste0(c("median_hh_", "gini_index", "pct_latino", 
                "pct_black", "pct_povert", "pct_emp_ed_health", "pct_emp_hospitality",
                "pct_emp_govt", "pct_bachel", "under18", "over65",
                "pop_densit", "ln_pop", "pritz18pct"), 
              collapse = " + ")

# Run models 
modelform <- formula(paste("tax_turnout ~ cases_1k + factor(COUNTYNAME) + ", cov, sep = ""))
mod1 <- lm_robust(modelform, data = df, se_type="HC2")
mod1b <- lm(modelform, data = df)

## Chicago deaths 
modelform <- formula(paste("tax_turnout ~ deaths_1k + ", cov, sep = ""))
mod2 <- lm_robust(modelform, data = df, se_type="HC2")
mod2b <- lm(modelform, data = df)

## Employment and earnings change 
modelform <- formula(paste("tax_turnout ~ emp_pct_loss + factor(COUNTYNAME) + ", cov, sep = ""))
mod3 <- lm_robust(modelform, data = df, se_type="HC2")
mod3b <- lm(modelform, data = df)

modelform <- formula(paste("tax_turnout ~ pay_pct_loss + factor(COUNTYNAME) + ", cov, sep = ""))
mod4 <- lm_robust(modelform, data = df, se_type="HC2")
mod4b <- lm(modelform, data = df)

## Lockdowns 
# (No FE b/c restrictions are constant within counties)
modelform <- formula(paste("tax_turnout ~ treated_lock + border_seg_fe + ", cov, sep = ""))
mod5 <- lm_robust(modelform, data = df, se_type="HC2")
mod5b <- lm(modelform, data = df)

# Produce table
stargazer::stargazer(mod1b, mod2b, mod3b, mod4b, mod5b,
                     se = c(starprep(mod1), starprep(mod2),
                            starprep(mod3), starprep(mod4),
                            starprep(mod5)), 
                        type = "text", 
                     keep = c("case", "death", "lock", "mean", "emp_pct", "pay_pct"),
                     star.cutoffs = 0.05,
                     omit.stat = "f")

## TABLE A7 --------------------------------------------------

cov <- paste0(c("median_hh_", "gini_index", "pct_latino", 
                "pct_black", "pct_povert", "pct_emp_ed_health", "pct_emp_hospitality",
                "pct_emp_govt", "pct_bachel", "under18", "over65",
                "pop_densit", "ln_pop", "pritz18pct"), 
              collapse = " + ")

# Run models 
modelform <- formula(paste("I(pritz_18_22_diff*100) ~ yestaxpct + factor(COUNTYNAME) + ", cov, sep = ""))
mod0 <- lm_robust(modelform, data = df, se_type="HC2")
mod0b <- lm(modelform, data = df)
# 1 p.p. increase in voting for tax reform -> 0.4 p.p. increase in Pritzker vote swing 

# Effects of pandemic on Pritzker vote shift
# Define covariate conditioning set
cov <- paste0(c("median_hh_", "gini_index", "pct_latino", 
                "pct_black", "pct_povert", "pct_emp_ed_health", "pct_emp_hospitality",
                "pct_emp_govt", "pct_bachel", "under18", "over65",
                "pop_densit", "ln_pop", "pritz18pct"), 
              collapse = " + ")

# Run models 
modelform <- formula(paste("I(pritz_18_22_diff*100) ~ cases_1k + factor(COUNTYNAME) + ", cov, sep = ""))
mod1 <- lm_robust(modelform, data = df, se_type="HC2")
mod1b <- lm(modelform, data = df)

# Deaths
# Define covariate conditioning set
cov <- paste0(c("median_hh_", "gini_index", "pct_latino", 
                "pct_black", "pct_povert", "pct_emp_ed_health", "pct_emp_hospitality",
                "pct_emp_govt", "pct_bachel", "under18", "over65",
                "pop_densit", "ln_pop", "pritz18pct"), 
              collapse = " + ")

modelform <- formula(paste("I(pritz_18_22_diff*100) ~ deaths_1k + ", cov, sep = ""))
mod2 <- lm_robust(modelform, data = df, se_type="HC2")
mod2b <- lm(modelform, data = df)

## Employment and earnings change 
modelform <- formula(paste("I(pritz_18_22_diff*100) ~ emp_pct_loss + factor(COUNTYNAME) + ", cov, sep = ""))
mod3 <- lm_robust(modelform, data = df, se_type="HC2")
mod3b <- lm(modelform, data = df)

modelform <- formula(paste("I(pritz_18_22_diff*100) ~ pay_pct_loss + factor(COUNTYNAME) + ", cov, sep = ""))
mod4 <- lm_robust(modelform, data = df, se_type="HC2")
mod4b <- lm(modelform, data = df)

## Lockdowns 
# (No FE b/c restrictions are constant within counties)
modelform <- formula(paste("I(pritz_18_22_diff*100) ~ treated_lock + border_seg_fe + ", cov, sep = ""))
mod5 <- lm_robust(modelform, data = df, se_type="HC2")
mod5b <- lm(modelform, data = df)

# Produce table
stargazer::stargazer(mod0b, mod1b, mod2b, mod3b, mod4b, mod5b,
                     se = c(starprep(mod0), starprep(mod1), starprep(mod2),
                            starprep(mod3), starprep(mod4),
                            starprep(mod5)), 
                       type = "text", 
                     keep = c("yes", "case", "death", "lock", "mean", "emp_pct", "pay_pct"),
                     star.cutoffs = 0.05,
                     omit.stat = "f")

## FIGURE A10 -------------------------------------------
# Source: https://www.covidstates.org/executive-approval
gov <- read.csv("./data/sm_data/Governor Approval.csv", na.strings = "") 
gov <- gov[,1:13] %>%
  mutate(val = coalesce(!!!select(., X20.Apr, X20.May, X20.Jun, X20.Jul, X20.Aug, 
                                  X20.Sep, X20.Oct, X20.Nov, X20.Nov.1)), 
         date = as.Date(End.Date, format = "%d-%b-%y")) %>%
  filter(Var == "Avg. Approval rating" & !is.na(val) & 
           date <= as.Date("2020-11-03")) %>%
  select(date, State, val) %>%
  mutate(val = as.numeric(gsub(pattern="%", x=val, "")), 
         State = ifelse(State=="National", "All states", State))

ggplot(gov) + 
  geom_point(aes(date, val/100, col=State)) + 
  # geom_label(aes(date, val/100, col=State, label=paste0(val, "%"))) + 
  geom_line(aes(date, val/100, col=State)) + 
  labs(x="", y="Percentage saying\nthey approve of\nthe way their\ngovernor is\nhandling the\nCOVID-19 outbreak", 
       col="") + 
  theme(axis.title.y = element_text(angle=0, vjust=0.5), 
        legend.position = ) + 
  scale_y_continuous(limits = c(0,.75), labels = scales::percent, 
                     breaks = c(0, 25, 50, 75)/100) +
  scale_color_manual(values=c("black", "red", "blue"))

# Alternative governor approval 
# Source: https://napequity.org/wp-content/uploads/NAPE_governor-approval-ratings_edits_5.21.20_vfinal-002.pdf
gov <- read.csv("./data/sm_data/gov_approval2.csv") %>%
  filter(!is.na(pre)) %>%
  mutate(diff = post-pre,
         diff2 = ifelse(diff>0, paste0("+",diff), diff)) %>%
  mutate(state = reorder(state, post)) %>%
  pivot_longer(cols = c(pre, post)) %>%
  mutate(name = factor(name, levels=c("pre", "post")))

ggplot(gov) + 
  geom_bar(aes(value, state, fill=name), 
           stat="identity", position = position_dodge(), 
           alpha=0.4) + 
  geom_bar(data=gov %>% mutate(value = ifelse(state %in% c("Illinois", "Arizona"), 
                                              value, NA),
                               diff = ifelse(state %in% c("Illinois", "Arizona"), 
                                             diff, NA)),
           aes(value, state, fill=name), 
           stat="identity", position = position_dodge()) +
  geom_text(data=gov %>% filter(name=="pre"),
            aes(88, state, label=diff2, col=diff)) + 
  scale_color_gradient(low="red", high="blue") + 
  scale_fill_manual(values=c("gray65", "gray20"),
                    labels=c("Pre-pandemic (Q4 2019)",
                             "Post-pandemic (Q1 2020)")) + 
  labs(y="", x="Governor approval\nrating", fill="") + 
  guides(col="none") + 
  theme(legend.position = "bottom")


## TABLE A8 -------------------------------------------

# Load in Arizona data 
df_az <- read.csv("./data/sm_data/master_zip_data_az.csv")

# Define covariate conditioning set
cov <- paste0(c("median_hh_", "gini_index", "pct_latino", 
                "pct_black", "pct_povert", "pct_emp_ed_health", "pct_emp_hospitality",
                "pct_emp_govt", "pct_bachel", "under18", "over65",
                "pop_densit", "ln_pop", "ducey18pct"), 
              collapse = " + ")

# Case counts
modelform <- formula(paste("yestaxpct ~ cases_1k + factor(countyname) + ", cov, sep = ""))
mod1 <- lm_robust(modelform, data = df_az, se_type="HC2")
mod1b <- lm(modelform, data = df_az)

# Effect of 1 s.d. increase
mod1$coefficients["cases_1k"] * sd(df$cases_1k, na.rm = TRUE)

## Employment and earnings change
modelform <- formula(paste("yestaxpct ~ emp_pct_loss + factor(countyname) + ", cov, sep = ""))
mod2 <- lm_robust(modelform, data = df_az, se_type="HC2")
mod2b <- lm(modelform, data = df_az)

modelform <- formula(paste("yestaxpct ~ pay_pct_loss + factor(countyname) + ", cov, sep = ""))
mod3 <- lm_robust(modelform, data = df_az, se_type="HC2")
mod3b <- lm(modelform, data = df_az)

# Produce table
stargazer::stargazer(mod1b, mod2b, mod3b,
                     se = c(starprep(mod1), starprep(mod2),
                            starprep(mod3)), 
                        type = "text", 
                     keep = c("case", "death", "lock", "mean", "emp_pct", "pay_pct"),
                     star.cutoffs = 0.05,
                     omit.stat = "f")

## TABLE A9 ----------------------------------------------
anes %>%
  filter(panel==1) %>%
  group_by(year) %>%
  summarise(yes_mil_tax = mean(yes_mil_tax, na.rm = TRUE),
            mil_tax = mean(mil_tax, na.rm = TRUE),
            working = mean(working, na.rm = TRUE),
            fin_worry = mean(fin_worry, na.rm = TRUE),
            covid_infect = mean(covid_infect, na.rm = TRUE),
            dem = mean(pid=="Democrat", na.rm = TRUE),
            ind = mean(pid=="Independent", na.rm = TRUE),
            rep = mean(pid=="Republican", na.rm = TRUE),
            egal_idx = mean(egal_idx, na.rm = TRUE)
  )

## TABLE A10 -------------------------------------

# Convert panel to first differences
anes_fd <- anes %>% 
  filter(panel==1) %>%
  # Rescale financial worry
  mutate(fin_worry_z = (fin_worry - mean(fin_worry, na.rm = TRUE)) / 
           sd(fin_worry, na.rm = TRUE)) %>%
  group_by(id, pid16) %>%
  summarise(mil_tax_fd = mil_tax[year==2020] - mil_tax[year==2016],
            yes_mil_tax_fd = yes_mil_tax[year==2020] - yes_mil_tax[year==2016],
            covid_infect = max(covid_infect),
            chg_working = working[year==2020] - working[year==2016],
            chg_worry =  fin_worry_z[year==2020] - fin_worry_z[year==2016]
  ) %>%
  ungroup()

# Estimate models 
mod1 <- lm_robust(mil_tax_fd ~ covid_infect,
                  data=anes_fd, se_type = "HC2")
mod1b <- lm(mil_tax_fd ~ covid_infect,
            data=anes_fd)

mod2 <- lm_robust(mil_tax_fd ~ chg_worry,
                  data=anes_fd, se_type = "HC2")
mod2b <- lm(mil_tax_fd ~ chg_worry,
            data=anes_fd)

mod3 <- lm_robust(mil_tax_fd ~ chg_working,
                  data=anes_fd, se_type = "HC2")
mod3b <- lm(mil_tax_fd ~ chg_working,
            data=anes_fd)

mod4 <- lm_robust(mil_tax_fd ~ covid_infect + chg_worry + chg_working,
                  data=anes_fd, se_type = "HC2")
mod4b <- lm(mil_tax_fd ~ covid_infect + chg_worry + chg_working,
            data=anes_fd)

stargazer::stargazer(mod1b, mod2b, mod3b, mod4b,
                     se = c(starprep(mod1), starprep(mod2),
                            starprep(mod3), starprep(mod4)), 
                     type = "text", 
                     omit.stat = "f")


## TABLE A11 ----------------------------
# Convert panel to first differences
anes_fd <- anes %>% 
  filter(panel==1) %>%
  # Rescale financial worry
  mutate(fin_worry_z = (fin_worry - mean(fin_worry[year==2016], na.rm = TRUE)) / 
           sd(fin_worry[year==2016], na.rm = TRUE),
         egal_idx_sd2016 = sd(egal_idx[year==2016], na.rm = TRUE), 
         egal1_sd2016 = sd(egal1[year==2016], na.rm = TRUE), 
         egal2_sd2016 = sd(egal2[year==2016], na.rm = TRUE), 
         egal3_sd2016 = sd(egal3[year==2016], na.rm = TRUE), 
         egal4_sd2016 = sd(egal4[year==2016], na.rm = TRUE) 
  ) %>%
  group_by(id, pid16) %>%
  summarise(mil_tax_fd = mil_tax[year==2020] - mil_tax[year==2016],
            yes_mil_tax_16 = yes_mil_tax[year==2016],
            yes_mil_tax_fd = yes_mil_tax[year==2020] - yes_mil_tax[year==2016],
            egal_idx_fd = (egal_idx[year==2020] - egal_idx[year==2016]),
            egal_idx_fd_z = (egal_idx[year==2020] - egal_idx[year==2016]) / egal_idx_sd2016,
            egal1_fd_z = (egal1[year==2020] - egal1[year==2016]) / egal1_sd2016,
            egal2_fd_z = (egal2[year==2020] - egal2[year==2016]) / egal2_sd2016,
            egal3_fd_z = (egal3[year==2020] - egal3[year==2016]) / egal3_sd2016,
            egal4_fd_z = (egal4[year==2020] - egal4[year==2016]) / egal4_sd2016,
            covid_infect = max(covid_infect),
            chg_working = working[year==2020] - working[year==2016],
            chg_worry =  fin_worry_z[year==2020] - fin_worry_z[year==2016]
  ) %>%
  ungroup() %>%
  distinct()

mod1 <- lm_robust(egal_idx_fd_z ~ covid_infect,
                  data=anes_fd, se_type = "HC2")
mod1b <- lm(egal_idx_fd_z ~ covid_infect,
            data=anes_fd)

mod2 <- lm_robust(egal_idx_fd_z ~ chg_worry,
                  data=anes_fd, se_type = "HC2")
mod2b <- lm(egal_idx_fd_z ~ chg_worry,
            data=anes_fd)

mod3 <- lm_robust(egal_idx_fd_z ~ chg_working,
                  data=anes_fd, se_type = "HC2")
mod3b <- lm(egal_idx_fd_z ~ chg_working,
            data=anes_fd)

mod4 <- lm_robust(egal_idx_fd_z ~ covid_infect + chg_worry + chg_working,
                  data=anes_fd, se_type = "HC2")
mod4b <- lm(egal_idx_fd_z ~ covid_infect + chg_worry + chg_working,
            data=anes_fd)

stargazer::stargazer(mod1b, mod2b, mod3b, mod4b,
                     se = c(starprep(mod1), starprep(mod2),
                            starprep(mod3), starprep(mod4)), 
                         type = "text", 
                     omit.stat = "f")

## TABLE A12 ----------------------------------------------

anes_fd <- anes %>% 
  filter(panel==1) %>%
  # Rescale financial worry
  mutate(fin_worry_z = (fin_worry - mean(fin_worry[year==2016], na.rm = TRUE)) / 
           sd(fin_worry[year==2016], na.rm = TRUE),
         egal_idx_sd2016 = sd(egal_idx[year==2016], na.rm = TRUE), 
         egal1_sd2016 = sd(egal1[year==2016], na.rm = TRUE), 
         egal2_sd2016 = sd(egal2[year==2016], na.rm = TRUE), 
         egal3_sd2016 = sd(egal3[year==2016], na.rm = TRUE), 
         egal4_sd2016 = sd(egal4[year==2016], na.rm = TRUE) 
  ) %>%
  group_by(id, pid16) %>%
  summarise(mil_tax_fd = mil_tax[year==2020] - mil_tax[year==2016],
            yes_mil_tax_16 = yes_mil_tax[year==2016],
            yes_mil_tax_fd = yes_mil_tax[year==2020] - yes_mil_tax[year==2016],
            egal_idx_fd = (egal_idx[year==2020] - egal_idx[year==2016]),
            egal_idx_fd_z = (egal_idx[year==2020] - egal_idx[year==2016]) / egal_idx_sd2016,
            egal1_fd_z = (egal1[year==2020] - egal1[year==2016]) / egal1_sd2016,
            egal2_fd_z = (egal2[year==2020] - egal2[year==2016]) / egal2_sd2016,
            egal3_fd_z = (egal3[year==2020] - egal3[year==2016]) / egal3_sd2016,
            egal4_fd_z = (egal4[year==2020] - egal4[year==2016]) / egal4_sd2016,
            covid_infect = max(covid_infect),
            chg_working = working[year==2020] - working[year==2016],
            chg_worry =  fin_worry_z[year==2020] - fin_worry_z[year==2016],
            health20 = health[year==2020] - 1,
            working20 = working[year==2020]
  ) %>%
  ungroup() %>%
  distinct()

mod1 <- lm_robust(yes_mil_tax_fd ~ covid_infect * working20,
                  data=anes_fd, se_type = "HC2")
mod1b <- lm(yes_mil_tax_fd ~ covid_infect * working20,
            data=anes_fd)

mod2 <- lm_robust(yes_mil_tax_fd ~ covid_infect * health20,
                  data=anes_fd, se_type = "HC2")
mod2b <- lm(yes_mil_tax_fd ~ covid_infect * health20,
            data=anes_fd)

stargazer::stargazer(mod1b, mod2b,
                     se = c(starprep(mod1), starprep(mod2)), 
                     type = "text", 
                     omit.stat = "f")


## FIGURE A11 -----------------------------------------

mod1 <- lm_robust(yes_mil_tax_fd ~ covid_infect * health20,
                  data=anes_fd, se_type = "HC2")
mod2 <- lm_robust(yes_mil_tax_fd ~ covid_infect * factor(health20),
                  data=anes_fd, se_type = "HC2")


margfx <- rbind(margins(mod1, variables = "covid_infect", 
                        at = list(health20=seq(0,4,length.out=20))) %>% 
                  summary() %>% 
                  mutate(type="linear"),
                margins(mod2, variables = "covid_infect", 
                        at = list(health20=0:4)) %>% 
                  summary() %>% 
                  mutate(type="categorical"))

ggplot(margfx) + 
  geom_hline(yintercept = 0, lty="dashed", col="grey") + 
  geom_point(data=margfx %>% 
               filter(type=="categorical"), 
             aes(health20, AME))  + 
  geom_linerange(data=margfx %>% 
               filter(type=="categorical"), 
             aes(health20, ymin=lower, ymax=upper))  + 
  geom_line(data=margfx %>% 
               filter(type=="linear"), 
             aes(health20, AME)) + 
  geom_ribbon(data=margfx %>% 
              filter(type=="linear"), 
            aes(x=health20, ymin=lower, ymax=upper), 
            alpha=0.3) + 
  scale_x_continuous(breaks=0:4, 
                     labels = c("Poor", "Fair", "Good", 
                                "Very good", "Excellent")) + 
  labs(x="Self-reported\nhealth in 2020", 
       y="Implied marginal\neffect of\nCovid-19\ninfection on\nsupport for\nprogressive\ntaxation") + 
  theme(axis.title.y = element_text(angle=0, vjust=0.5), 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))



